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ABSTRACT 

Any  change  in  neutronic  properties  in  a  reactor  operat- 
ing at  steady  state  will  result  in  a  change  in  the  equilib- 
rium neutron  flux  and  hence,  the  power  of  the  reactor.   A 
main  cause  for  a  change  in  neutronic  properties  is  the  high 
temperature  attained  in  a  reactor,  which  produces  a  feedback 
in  the  reactor  operation.   The  response  of  the  reactor  to 
a  particular  feedback  is  analyzed  by  using  Liapunov's 
Second  Method  to  specify  stability  regimes.   Both,  the 
point-kinetics  model  and  the  distributed  parameters  system 
are  analyzed . 

Data  from  a  typical  thermal  and  fast  reactors  is  used 
to  specifically  determine  the  stability  domains. 
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I .   LIAPUNOV'S  SECOND  METHOD 

A.   HISTORY 

In  1892  the  Russian  mathematician  A.  M.  Liapunov  pub- 
lished a  long  paper  dealing  with  the  problem  of  stability 
of  motion  (1_)  .   This  paper  was  translated  into  French  in 
1907  and  reprinted  in  America  in  1949.   Liapunov's  Theory 
received  by  then  only  little  attention  due  to  the  diffi- 
culty in  understanding  the  advanced  mathematical  theorems, 
to  the  abstract  way  it  was  presented  and  to  its  lack  of 
practical  application  and  for  a  long  time  it  was  nearly 
forgotten.   About  35  years  ago,  Soviet  mathematicians  re- 
sumed the  investigation  of  Liapunov's  Theory  and  its  excel- 
lent application  in  several  technical  fields,  mainly  in 
control  engineering,  was  noticed.   Significant  work  in  this 
area  was  published  by  Malkin  (_2)  ,  Letov  (3_)  ,  Lur 'e  (4)  and 
Chetaev  (_5)  .   The  excellent  paper  by  Massera  (6^)  and  the 
translation  to  English  of  most  of  the  Russian  works  stirred 
up  considerable  curiosity  in  this  country.   Bellman  (7_)  in 
1953  published  an  excellent  work  concerning  stability,  but 
the  section  on  Liapunov's  method  is  difficult  to  comprehend 
According  to  this  author,  Hahn  (J3)  ,  Krasovskii  (9_)  and 
LaSalle  and  Lefschetz  (10.)  have  the  best  treatises  on  the 
subject,  the  consolidation  of  the  concepts  of  Liapunov's 
Theory  and  the  clear  presentation  of  it,  make  these  three 
books  the  main  references  for  workers  in  the  field. 
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Recently,  Zubov.  (1_1)  found  the  best  extension  of  the 
Liapunov's  Theory,  to  include  the  analysis  of  partial  dif- 
ferential equations,  which  is  by  now  the  main  topic  for 
research.   Yet,  a  good  amount  of  work  remains  to  be  done. 
All  the  previous  works,  including  Liapunov's  original 
theory,  were  devoted  to  the  analysis  of  stability  of 
ordinary  differential  equations. 

Many  authors  have  applied  Zubov's  extension  to  works 
concerning  vibrations,  reactor  physics,  hydrodynamics, 
magne tohydr odynamics  and  control  processes.   The  list  of 
references  will  provide  the  interested  reader  with  the 
main  works  in  this  field. 

B.   THEORY 

The  name  "second  method"  (or  "direct  method")  is  of 
historical  origin.   Liapunov's  Theory  is  divided  in  two 
categories.   The  "first  method"  which  comprises  all  pro- 
cedures in  which  the  explicit  form  of  the  solutions  is 
used,  specially  when  represented  by  infinite  series  and  the 
investigation  of  stability  "in  the  small"  (local  stability) 
studies  the  singular  points  of  a  nonlinear  differential 
equation  by  using  the  appropriate  linearized  version  of  the 
differential  equation  near  the  singular  point.   The  "second 
method"  attempts  to  make  stability  statements  by  using  (in 
addition  to  the  differential  equations)  suitable  functions, 
called  Liapunov  functions,  which  are  defined  in  the  motion 
space.   This  method  which  deserves  special  study  due  to 
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its  inherently  advantages  over  most  of  the  conventional 
methods,  investigates  the  stability  "in  the  large"  (global 
stability) . 

Liapunov's  Second  Method  is  in  essence  a  more  general 
expression  of  the  Hamiltonian  (total  energy),  and  it  is 
based  in  the  statement  that  a  physical  system  loses  poten- 
tial energy  in  a  neighborhood  of  a  point  of  stable  equilib- 
rium.  It  is  said  that  it  is  more  general  than  the  total 
energy  method,  because  unlike  the  energy  of  a  system,  the 
Liapunov  function,  denoted  V,  is  not  unique.   This  is  the 
main  reason  why  the  "second  method"  is  a  tool  in  the  analy- 
sis of  stability  of  dynamical  systems,  which  must  be  used 
with  considerable  skill.   One  of  the  main  features  of  this 
method  is  its  appeal  to  geometric  intuition;  V(x,t)  can  be 
seen  as  a  measure  of  the  "distance"  of  the  state  (x,t)  from 
the  origin,  in  the  state  space  and  the  variation  of  this 
"distance"  (norm  of  the  differential  equation)  as  t  varies 
will  provide  definite  bounds  of  stability  regions  for  the 
prescribed  system  under  consideration. 

Normally  stability  with  respect  to  one  norm  does  not 
imply  stability  with  respect  to  another  norm.   This  diffi- 
culty does  not  appear  in  finite-dimensional  systems  since 
all  norms  defined  on  a  finite-dimensional  vector  space  are 
equivalent . 

1 .   Fundamental  Concepts  and  Definitions 

Stability  is  a  property  of  certain  systems  of  dif- 
ferential equations  and  is  basically  concerned  with  the 
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question  of  whether  or  not  the  dynamic  system  will  return 
to  a  particular  state  after  it  has  been  disturbed  in  some 
way  . 

Differential  equations  in  their  most  general  way 
can  be  expressed  as 

iaiJLti  .  F[u(-jt)it] 

where  U(x,t)  is  a  multidimensional  function  of  space  and 
time.   It  may  well  happen  that  F  depends  upon  U(t)  alone 
and  not  upon  time  explicitly.   Then  the  previous  equation 
assumes  the  simpler  form 

U(t)  =  F[U(t)] 
A  system  of  this  nature  is  known  as  autonomous.   Since  it 
is  not  intended  in  this  work  to  expose  the  reader  with  the 
mathematical  aspect  of  Liapunov's  Theory,  in  general,  the 
presentation  of  the  theory  will  be  made  without  proofs. 
Let  us  denote  f2(R)  the  spherical  region  ||x||  <  R  and  by 
A(R)  the  sphere  boundary  ||x|[  =  R.   The  matter  of  concern 
is  the  stability  of  the  origin.   Initiating  the  motion  at 
a  point  x  ,  the  origin  is  said  to  be: 

(a)  Stable  whenever  the  path  remains  in  the  spherical 
region  fi(R);  that  is, the  path  never  reaches  the 
boundary  sphere  A(R),  Figure  1. 

(b)  Asymptotically  stable  whenever  it  is  stable  and 
in  addition  the  path  tends  to  the  origin  as  time 
increases  indefinitely. 
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Unstab le 


A(R) 


Figure  1.   Domains  of  stability. 
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(c)   Unstable  whenever  the  path  reaches  the  boundary 
sphere  A (R) . 
For  autonomous  systems: 

A  Liapunov  function,  V(U),  is  defined  as  a  posi- 
tive definite  scalar  function  with  the  following  properties 

(a)  V(U)  is  continuous  with  its  first  partial  deriva- 
tives in  a  certain  open  region  about  the  origin 

(b)  V(o)  -  0  only 

(c)  Outside  the  origin  and  always  inside  the  open 
region,  V  is  non-negative  and  vanishes  only  at  the 
origin.   The  origin  is  an  isolated  minimum  of  V; 
and  in  addition 

V  <  0  in  the  open  region 

If  V  =  0,  the  stability  is  neutral 

If  V  <  0,  the  stability  is  asymptotic. 

For  nonautonomous  systems: 

Let  us  define  V. (U)  as  previously.   A  Liapunov 
function,  V(U,t)  is  defined  as  a  positive  definite  scalar 
function  with  the  following  properties: 

(a)  V(U,t)  is  defined  in  the  open  region  for  all 
t  >  0. 

(b)  V(o,t)  =  0  for  t  >  0. 

(c)  V1(U)  <^  V(U,t)  for  all  x  in  the  open  region  and 
all  t  >^  0  ;  and  in  addition 

V  <  0  in  the  open  region 
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where 


Summarizing  these  conditions,  the  three  main 
steps  to  test  if  a  scalar  function  is  a  Liapunov  func- 
tion a  r  e  : 

1.  V  >  0 

2.  V(o)  =  0 

3.  V  <  0 
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2 .   Methods  to  Construct  a  Liapunov  Function 

Several  authors  have  attempted  to  present  guide- 
lines for  generating  Liapunov  Functions.   Unfortunately, 
most  of  these  methods  are  extremely  restrictive,  limita- 
tions being  imposed  primarily  by  the  number  of  nonlinear 
terms  and  the  order  of  the  system.   Ingwerson   (JJO  pre- 
sents a  Table  to  generate  Liapunov  Functions  for  linear 
autonomous  differential  equations  up  to  the  fourth  order. 
Barbashin  (1_3)  and  Lur'e  and  Rozenvasser  (1J0  tried  to 
establish  some  rules  for  construction  of  Liapunov  Functions 
Schultz  (15 )  has  one  of  the  best  methods  available  now;  it 
is  called  the  Variable  Gradient  Method,  being  the  least 
restrictive  for  the  case  of  autonomous  systems.   Its  only 
restriction  is  that  all  nonlinear i ties  must  be  single- 
valued  . 

In  general,  literature  late  in  1961  and  1962  became 
saturated  with  methods  for  generating  Liapunov  Functions, 
most  of  them  with  the  restrictions  previously  mentioned. 
Furthermore,  such  methods,  when  applicable,  were  directed 
to  autonomous  systems  only. 

Although  development  of  Liapunov  stability  theory 
and  applications  to  ordinary  differential  equations  has 
progressed  rapidly,  its  application  to  partial  differential 
equations  has  remained  limited.   The  importance  of  partial 
differential  equations  in  the  fields  of  reactor  physics, 
hydrodynamics,  control  processes,  etc.  has  motivated  inves- 
tigations of  possible  ways  to  extend  Liapunov  stability 
theory  to  partial  differential  equations. 
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In  this  study  Functional  Analysis  plays  an  important 
role,  because  the  stability  of  the  equilibrium  solution  is 
defined  in  terms  of  the  norm  induced  by  the  inner  product  of 
the  Hilbert  space  on  which  the  solutions  of  the  system  are 
defined.   Thus,  with  proper  choice  of  inner  product  or  norm, 
the  square  of  the  norm  becomes  the  Liapunov  functional  which 
establishes  asymptotic  stability.   In  fact  what  is  being 
done  is  the  construction  of  a  scalar  function  of  the  distance 
between  the  solution  and  the  equilibrium  point  of  interest. 
If  this  distance  function,  evaluated  along  the  solution 
paths,  obeys  certain  inequalities,  then  statements  concern- 
ing the  stability  of  the  equilibrium  point  can  be  made. 
Movchan  (1_£)  defines  stability  in  terms  of  two  metrics, 
rather  than  one,  to  be  more  restrictive  on  the  initial 
states.   The  choice  of  the  initial  state  space  and  the 
metric  is  crucial  in  the  formulation  of  stability  problems 
in  the  framework  of  Liapunov's  Second  Method.   For  cases  in 
which  only  the  behavior  of  some  of  the  state  variables  is 
of  importance  or  some  function  of  the  state  variables  is  of 
interest  to  the  analyst,  then,  for  these  cases,  it  is  mean- 
ingful to  define  Liapunov  stability  with  respect  to  two 
metrics.   For  instance,  one  may  be  only  interested  in  the 
behavior  of  the  maximum  deflections  in  an  elastic  system, 
regardless  of  the  velocity  and  acceleration  involved  into 
attained  it . 

According  to  Kastenberg  &  Ziskind(l_7^)  extreme  care 
has  to  be  taken  when  interpreting  the  obtained  results, 
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because  severe  peaking  in  some  of  the  state  variables  of  the 
system  can  cause  the  L~  norm  of  the  system  trajectory  to 
move  out  of  the  domain.   Wang  (18)  and  Parks  (19_)  have  de- 
voted special  attention  to  the  study  of  panel  flutter  which 
represents  a  linear  distributed  parameter  system.   It  is 
concluded  that  a  good  comprehension  of  Hilbert  and  Banach 
spaces  will  provide  the  analyst  with  an  excellent  tool  for 
generating  Liapunov  functionals. 

3 .   Comparison  with  Other  Methods 

The  advantages  of  Liapunov's  Second  Method  over  the 
conventional  methods  used  for  stability  analysis  can  be 
summarized  as  follows: 

(a)  It  employs  the  system  equations  without  resorting 
to  approximations.   Normally  a  distributed  system  is  approxi- 
mated by  a  lumped  parameter  model  having  finite  or  infinite 
number  of  degrees  of  freedom.   This  method  is  not  satisfac- 
tory because  quite  often  it  exhibits  characteristics  which 
do  not  agree  with  the  physical  nature  of  the  problem. 
Liapunov's  method  is  in  general  a  more  reliable  method. 

(b)  There  is  no  theoretical  limit  on  the  number  of  non- 
linearities  or  on  the  order  of  the  differential  equation  to 
which  Liapunov's  Second  Method  can  be  applied. 

(c)  The  laborious  work  implied  in  finding  the  system 
solution  is  avoided.   Normally  the  system  equations  are  in- 
tegrated numerically  for  some  given  perturbations  in  the 
initial  conditions.   This  method  is  time  consuming  and 
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sometimes  does  not  give  an  accurate  presentation  of  the 
system  behavior. 

(d)  The  relationship  between  the  system  stability  and 
the  system  parameters  is  directly  extracted  from  the  analy- 
sis of  Liapunov's  method  conditions  for  stability. 

(e)  There  are  no  mathematical  restrictions  with  res- 
pect to  uniqueness. 

The  method  presents  also  some  deficiencies;  among 
them  the  main  ones  are: 

(a)  The  construction  of  the  required  Liapunov  Function 
for  the  system  under  consideration  is  without  doubt  the 
most  difficult  part  of  the  task  due  to  the  fact  that  there 
is  no  guideline  to  construct  it  and  success  depends  mainly 
upon  the  ingenuity  and  experience  of  the  analyst. 

(b)  The  interpretation  of  the  obtained  results  has  to 
be  done  carefully. 

(c)  Integral  inequalities  play  an  important  role  in 
deriving  conditions  for  stability.   As  a  consequence,  the 
regions  of  stability  may  be  somewhat  looser  than  those  ob- 
tained by  other  methods. 

In  conclusion,  in  spite  of  the  difficulties  cited 
above,  the  advantages  of  the  Liapunov's  Second  Method  makes 
it  an  excellent  method  for  studying  stability. 

C.   APPLICATION  TO  NUCLEAR  REACTOR  CONTROL 

Liapunov's  Second  Method  has  been  applied  with  success 
by  Hsu  (20)  who  considers  the  linear  and  non-linear  cases, 
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analyzing  them  by  means  of  the  Spectral  Analysis  and  Liapu- 
nov's  Method,  and  comparing  the  results  obtained  by  each 
method.   llsu  shows  that  the  Spectral  Analysis  provides  a 
little  more  information  than  Liapunov's  Method,  but  this 
latter  eliminates  the  calculation  of  the  eigenvalues  of  the 
system,  and  for  the  nonlinear  system  analysis  bounds  of 
stability  are  presented. 

Kastenberg  (2JL)  presents  a  clear  comparison  of  Liapunov's 
Second  Method  with  the  Semigroup  Method,  and  the  Comparison 
Function  and  Maximum  Principle  Techniques.   In  the  Liapunov 
and  Semigroup  methods,  one  must  obtain  an  a  priori  bound  on 
the  system  nonlinearity  and  then  proceed.   When  employing 
the  maximum  principle  or  the  comparison  function  technique, 
an  a  priori  bound  on  the  system  nonlinearity  is  not  always 
required.   In  contrast,  one  must  give  a  bound  on  the  initial 
condition.   For  cases  which  are  stable  in  a  global  sense, 
the  results  of  the  various  methods  coincide. 

The  application  of  Liapunov's  Second  Method  to  nuclear 
reactor  control  seems  to  be  excellent,  mainly  for  the  spa- 
tially-dependent reactor  system,  due  to  the  fact  that  the 
method  allows  the  inclusion  of  any  number  of  nonlinear ities . 
This  will  permit  the  study  of  spatial  effects  of  temperature, 
the  main  feedback  effect  in  a  reactor,  control  rod  motion 
and  several  other  processes  within  a  reactor.   Also,  the 
inclusion  in  the  analysis  of  the  different  reactivity  coef- 
ficients, as  doppler  resonance  broadening  and  structural 
expansion  is  allowed. 
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As  pointed  before,  the  interpretation  of  the  results 
must  be  done  very  carefully  due  to  the  fact  that  a  Liapunov 
Function  is  not  unique. 
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II.   MATHEMATICAL  MODEL  OF  THE  REACTOR  SYSTEM 

Both  Thermal  and  Fast  reactors  are  described  essential- 
ly by  the  same  basic  dynamic  principles,  regardless  of 
material  and  geometry  considerations. 

For  a  given  reactor  size,  the  reactor  reactivity  depends 
on  the  neutron  cross-sections  and  on  the  relative  amounts 
and  densities  of  different  materials.   All  of  these  being 
affected  by  the  temperature,  the  reactivity  will  be  strong- 
ly coupled  to  the  power  of  the  reactor  and  the  reactor 
governing  equations  become  nonlinear. 

The  Thermal  Reactor  stability  is  analyzed  by  means  of 
the  lumped-parameter  (point-kinetics)  model,  in  which  the 
partial  differential  equations  are  reduced  to  ordinary  dif- 
ferential equations  via  spatial  discretization.   This  is 
justified  only  when  infinitesimal  small  perturbations  are 
introduced,  as  in  the  case  when  k  is  very  near  to  unity 
(very  small  departures  from  "critical").   Thermal  reactors 
are  generally  smaller  in  size  than  fast-breeder  reactors. 
The  Fast-Breeder  Reactor  stability  is  analyzed  using  the 
governing  partial  differential  equations  without  resorting 
to  any  type  of  approximations.   Due  to  possible  stronger 
space  effects  in  fast  reactors  than  in  thermal  reactors, 
the  space  and  time  of  the  state  variables  of  the  system 
must  be  maintained  during  the  analysis.   The  reactivity  in- 
sertion in  the  system,  Figure  2,  can  be  either  positive  or 
negative.   Both  cases  will  be  considered  in  this  work. 
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A.   POINT-KINETICS  MODEL 

1 .   Diffusion  Equation 

According  to  Meghreblian  and  Holmes  (2_2_)  ,  the  neu- 
tron diffusion  equation  is  written  in  the  following  way: 


D    V2f(r\t)    -    Ea*(r,t)    =   I  !A£i±l  - 

-    (l-3)v    Zf    pg    <J>(r,t)    -    Pg^Ai    C±(r.t) 


(1) 


Expressing  the  equation  for  a  slab  reactor  and  using  only 
one-delayed  neutron  group: 


D 


924>(x,t) 

a  2 
dx 


1  3<j)(xtt) 


MC*.t)  -  7   3T 


-  (l-3)v  Zf  pg  <J)(x,t)  -  pg  A  C(x,t) 


(2) 


It  is  looked  for  a  solution  which  is  separable  in  time  and 

space : 

<|>(x,t)  -  *(t)  F(x)  (a)   (3) 

C(x,t)  =  C(t)  F(x)  (b) 

Applying  relations  (3)  to  Equation  (2)  yields 


D  9    F(x) 

F(X>       8x2 


-    E 


4>(t) 


£  i(t)    - 


-  (i-e)v  zf  pg  <j>(t)  -  pg  a  c(t) 


Thus 


] 


1  92F(x)  1       )  _ 

F(x)  .2  D       /     a         4> 

dx  » 


-    (l-B)v    E,    pg    (J)(t)    -    pg    A    C(t) 


TO       [v   * 

]| 


4><t)    - 


-    -    B 


(4) 


(5) 
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From  Equation  (5)  it  is  obtained  that: 
32F(x)  _,   2 


+  B  F(x)  =  0 


(a) 


3x' 


(6) 


(DB2  +  Za)cJ)(t)  =  -  ^i(t)  +  (l-3)vEfPg<J)(t)+pg  X    C(t)    (b) 


Using  the  following  r elat ionships : 
,2     D 


k  = 


V  Ef  pg 


2  2 
E  (1+L  B  ) 
a 


£  = 


v(E  +DB  ) 
a 


and  k  =  Ak  +  1,  Equation  (6b)  becomes 


ti(t)-  [(l-B)Ak-3]4.(t)  +   ps  o  X  c(t)  (7) 

E  +DB 
a 

Let  P(t)  be  the  time  behavior  of  the  reactor  power 

generation 

p(t)  -  e  Z£  cj)(t)  (8) 


then    P(t)  =  £  E,  (f)(t),  and  the  equation  describing  the 


reactor  power  is 


E  Ef  Pg 

£P(t)=  [(l-3)Ak-3]P(t)  + J  X    C^) 

E  +  DB 
a 


(9) 


It  will  be  seen  in  the  next  paragraph  that  for  certain 
kinetics  problems  a  simplification  can  be  obtained  by 
assuming  the  infinite  delay  time  in  delayed  neutron  emis- 


sion : 


AC(o)  =  ~£  P(o) 
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This  expression  is  used  as  a  first  estimate  for  the  precur1 

sor  concentration. 

Then 


£P(t)  =  [Q-3)Ak-3]  P(t)  +  [1+Ak]  3  P(o) 


(10) 


2 .   Concentration  of  Precursors  Equation 

1^  C.(r,t)  =  -  X.  C±(r,t)  +  v3±  Ef  (j»(r,t)  (11) 

For  a  slab  reactor  and  using  one-delayed  neutron  group, 
Equation  (11)  can  be  expressed  as 

|^  C(x,t)  =  -  X  C(x,t)  +  v  3  2f  4>  (x,t)  (12) 

Using  Equations  (3)  to  separate  variables,  Equation  (12) 
becomes 

C(t)  =  -  X  C(t)  +  v  3  Ef  (f)(t)  (13) 

Then,  using  Equation  (8)  yields 


C(t)  =  -  X  C(t)  +  |  3  P(t) 


(14) 


From  Equation  (1A),  the  steady  state  concentration  of  the 


precursors  can  be  obtained  as: 

XC(o)  =  ^  P(o) 


(15) 


This  result  was  previously  used. 
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3  .   Energy  Equation 

Expressing  H(r,t)  as  the  energy  content  per  unit 
volume,  then  the  time  rate  of  change  of  H  must  be  given 
by  the  net  energy  gain  per  unit  time  and  volume  from  the 
fission  reactions  and  the  reactor  cooling: 

H(r,t)  =  pC  [T(r,t)  -  TQ(r)]  =  pCp  9(r,t)  (16) 

and 


|y  H(r,t)  -  £  Zf  <f>(r,t)  "  (?(.*, t) 
For  a  slab  reactor,  the  energy  balance  becomes 


(17) 


pCp  ~  6(x,t)  =  e  Ef  6(x,t)  -^(x,t) 

Let         6(x,t)  =  9(t)  F(x) 

<Kx,t)  =  <KO  F(x) 

(?  (x,t)  =(?(t)  F(x) 


Then 


pCp  6(t)  -  E  Z£  (})(t)  -    (t) 


(18) 

(a)   (19) 

(b) 

(c) 

(20) 


Let  Cf  =  pC   and  using  Equation  (8),  Equation  (20)  becomes 


C  0(t)  =  P(t)  -  (?(t) 


(21) 


Three  cases  are  normally  encountered  for  the  power  removal: 

(a)  Adiabatic,  (P   (t)  =0 

(b)  Constant  power  removal,  (p   (t)  =  A 

(c)  Newton's  law  of  cooling,  (p  (t)     =  hG(t) 
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4.   Feedback  Models 


a.   Temperature  Feedback  Model  //l 

The  multiplication  factor  is  temperature  and 
time  dependent 

k  =  k(T)  =  k[T(t)  ] 


Expanding  it  in  powers  of  [T(t)  -  T  ],  yields 


k(T)  = 


oo     n 

dk(T) 


z 


n  =  0 


[T(t)  -  T  ] 
d  t  o 


11 


thus 


n  =  0 


k(T)  =  kQ  [1  +  X  A  An  (T  -  Tq)] 


which  becomes,  upon  dropping  all  terms  beyond  n  =  1, 

k  =  k  [1  +  y(T-T  )]  =  (1  +  Ak  )(1  +  y0) 
or 


Ak(t)  =  Ak   +  yG(t) 


(22) 


(23) 


(24) 


(25) 


where  Ak   is  the  change  in  k  applied  to  the  reactor  at 
t  =  0  and  y  represents  the  temperature  coefficient  of 
reactivity,  which  could  be  either  (+)  or  (-) . 
b.   Feedback  Model  #2 

One  of  the  major  considerations  to  be  made  when 
accounting  for  the  temperature  dependence  of  the  multipli- 
cation is  that  described  by  the  nuclear  Doppler  effect. 

Thompson  and  Beckerley  (23)  expresses  that 
dependence  as 


dk  _    ,300, n 
dT    YD^  T  ; 


(26) 
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where  T  is  in  °K. 

From  (26),  Ak  can  be  expressed  explicitly  as 


Ak  -  Ak   +  ^-[K^V  -  I  <^V] 
o    1-n     t        o   1 0 


(27) 


being  n  =  1  a  special  case,  for  which  a  logarithmic  expres- 
sion is  found. 

Typical  values  of  n  are  1/2,  1  and  3/2. 


B.   DISTRIBUTED-PARAMETER  REACTOR  SYSTEM 
1 .   Diffusion  Equation 

The  diffusion  equation  is  used  in  this  case  without 
resorting  to  any  approximation,  then 


DV2cKr,t)  -  E  <Kr,t)  =  ~   |£  (r,t)  - 
a  vdt 

-(l-g)vZfPg      (j)(r,t)  -  pg  £X±  C.(r,t) 


(28) 


For  a  slab  reactor  and  using  one-group  delayed  neutron,  the 

diffusion  equation  takes  the  following  form: 

D^|(x,t)  -  Za  <Kx,t)  =  ~   f£<*-0  " 
8x 


-(l-3)v  Ef  pg  (J)(x,t)  -  pg  A  C(x,t) 


(29) 


The  cross  sections  should  be  written  in  a  proper  way  as 

time-dependent  but  these  variations  are.  assumed  negligible 

in  comparison  to  the  rapid  transients  in  $ ,  C  and  T.   The 

equation  for  a  non-trivial,  steady  state  solution  4>  (x)  , 

C  (x)  is: 
o 


92<l>  (x) 


D- 


9x' 


-  Z J    (x)  =  -(1-B)v£,  pg  <j>  (x)-pg  X    C  (x) 

3.   O  I  O  O 


(30) 
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Now  subtracting  Equation  (30)  from  Equation  (29)  yields 


D^|(x,t)  -  la    l|><x,t)  -  J  |fU.t)  - 

3x 

-(l-3)v  Z   pg  ip(x,t)  -  pg  X  ^  (x,t) 


where 


^(x,t)  =  f(x,t)  -  4>0(x) 
f (x,t)  =  C(x,t)  -  C  (x) 


(31) 

(a)   (32) 
(b) 


This  change  of  variables  moves  the  equilibrium  state  of  the 


system  from  (<i>  ,  C  )  to  (0,  0),  Figure  3. 
o    o 


Define 


and  using 


k  = 


V  Ef  pg 
a   1 


2   2 
a   =  1  +  L   B 


(33) 
(34) 


the  diffusion  equation  becomes 


Di_li2Lill  _  £   ij>(x,t)  +  (l-3)k  Z   a  i^(x,t)  + 


3x' 


^  ■/"    /        ^  \     1  3^(x,t) 

+  pg  X   if>   (x,t)  =  -  *dt>    / 


(35) 


Dimensionless  space  and  time  expressions  are  introduced  at 
this  point : 


y  = 


T  =  vE  t 
a 


(x  >  >  L) 


(a)   (36) 
(b) 


Then  Equation  (35)  becomes 


B2<Ky>T)    ,     r/,    ON1       ,  n  ,  /        s  ,    ka        ,   a ,        .        3^(y,T) 
— rw  +    [  (l-B)ka-l]^(y,T)+  ^        *p(y>"0    =      r^? 


8y 


(37) 


The  boundary  conditions  are:  (a)  ^(W,t)  =  0  (38) 

(b)i/;(-W,T)  =  o 
x 
where  W  =  -r—  expresses  the  dimensionless  half-width  of  the 
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Reactivity 
insert  ion 


T-   C 


(o,0) 


Figure  3.   Translation  of  equilibrium  states 
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slab  reactor  including  the  extrapolation  distance.   The 
initial  conditions  are:    (a)  \J;(y,o)  =  0  (39) 

(b)-^(y,o)  -  0 
2 .   Concentration  of  precursors  equation 


3C. (r , t) 

Tt 


=  -AC.  (r,  t)  +  v3.Ef  <J>(r,  t) 


(40) 


For  a  slab  reactor  and  using  one-group  delayed  neutron, 
Equation  (40)  is  written  in  the  following  way: 


5C(x, t) 

dt 


\C.(x,  t)  +  v3  £f  (J>(x,  t) 


(41) 


The  same  assumption  used  in  the  diffusion  equation,  with 
respect  to  the  cross  sections  is  used  here.   A  non  trivial, 


steady  state  solution  (j)  (x),  C  (x)  is  assumed: 

0  =  -  AC  (x)  +  V6  Z_  (J)  (x) 
o  to 

Subtracting  Equation  (42)  from  Equation  (41)  yields 

IT^LlLL   =  -x$   (x,t)  +  v3  Sf  *<x.t) 
Using  Equations  (36),  Equation  (43)  becomes 


(42) 


(43) 


^ii^ii.  -  ^f(y.T) 


+  vB  S.  ip<y,T) 


(44) 


The  initial  condition  for  the  equation  is 

^  (y,o)  =  0  (45) 

Introducing  the  multiplication  factor  in  Equation  (44)  yields 


vZ 


\±  (y.T) 

8t 


3aE  k 
-  Xf  (y,T)  +  —f-   *(y.T) 

t  o 


(46) 
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X 


Figure  A.   Dimensional 
slab  reactor 


geome  tr ical 
boundary 


^""  extrapolated 
di  stance 


Dimensionless 
slab  reactor . 
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3  .   Energy  Equation 

The  energy  balance  condition  is  normally  expressed 


as  : 


Energy  produced  -  Energy  removed  = 


Energy  stored  in 
the  system 


The  energy  produced  is  due  to  fission:   e  Zf  <j>(r,t) 

The  energy  removed  is  due  to  the  coolant:  (p    (r,t) 

Three  cases  are  normally  encountered  for  the  energy  removed 

from  a  system: 

(a)  Adiabatic  case,  (r    (r,t)  =  0 

(b)  Constant  power  removal,  for  which  the  steady 
state  requirements  are  normally  used,  (j     (r,o)=  eZ  (f>(r,o) 

or    (p   (r)  -  cEf0o(r) 

(c)  Newton's  law  of  cooling,  (p   (r,t)  =  h  0(r,t) 
It  is  noticed  that  these  cases  are  for  stationary-fuel 
reactors.   The  energy  equation  can  be  written  as: 


Let  C 


3Tf (r,t) 


8t 


E  £_  (j)(r,t)  -  (p   (r,t) 


(47) 


p,.C    be  the  heat  capacity  of  the  reactor  system 


expressed  as  energy  per  unit  volume  per  unit  temperature. 
The  subscripts  f  and  C  stand  for  fuel  and  coolant  respec- 
tively.  Using  the  same  assumption  as  used  previously  with 
respect  to  the  cross  section,  and  expressing  Equation  (47) 
for  a  slab  reactor,  yields 


9T  (x,t) 
Cf  ~ =  £  Zf  <J»(x,t)  -  (P   (x,t) 


(48) 
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a.   Adiabatic  Case 


3T  (x,t) 
Cf  — ^ -  e  Sf  *(x,t) 


(49) 


This  case  will  not  be  considered  in  this  work  because  it 
represents  accident  conditions. 

b.   Constant  Power  Removal 


Cf  -|^ =  E  Zf  Mx,t)  -    <P 


(50) 


Assuming  a  steady-state,  nontrivial  solution  cj>  (x)  it 
becomes 

0  =  e  If  <J>(x,t)  -  CP 

Subtracting  Equation  (51)  from  Equation  (50)  yields 


(51) 


m^i  __  ,      tu§t) 


(52) 


where        6(x,t)  -  Tr(x,t)  -  T.  (x) 

1  to 

Introducing  Equations  (36),  Equation  (52)  becomes 


(53) 


CfvZa  97  ^y'T)  =  c  Zf  My.T) 


(54) 


with  initial  condition: 

6(y,o)  =  0  (55) 

c.   Newton's  Law  of  Cooling 

3T  (x,t) 
Cf  ~ =  E  Ef  <J)(x,t)  -  h(Tf-  Tc)        (56) 

Assuming  a  steady-state,  nontrivial  solution,  <j)  (x)  and 
T_  (x) ,  Equation  (56)  becomes 


Cf  36^>t:)  -  E  Sf  *(x,t)  -  h6(x,t) 


(57) 
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Note  that  T   and  h  remain  constants.   This  is  explained  by 
c 

means  of  the  one  effective  coolant  temperature  assumption. 
Introducing  the  dimensionless  expressions  for  space  and  time, 
Equation  (62)  becomes: 


cf  v  za  363t,T)  =  c  zf  ^(y»T)  -  hQ(y>T) 


with  initial  condition:   6  (  y  >  o  )  =  0 
A.   Feedback  Models 


(58) 


(59) 


a.   Temperature  Feedback  Model  //l 

Through  the  temperature,  the  multiplication 
factor  is  now  dependent  on  both  space  and  time. 

k  =  k(T)  -  k[T(x,t) ]  (60) 

Expanding  it  in  powers  of  [T(x,t)  -  T  (x)]  yields: 


oo     n 

k(T)  =£   3  Up     [T(x,t)  -  T  (x)]n 
n  =  0 


Thus 

oo 

k(T)  -  k  [1  +Y       A  (T  -  T  )n] 
v/      o     ^    n        o 

n  =  l 

which  upon  dropping  all  terms  beyond  n  =  l,  reduces  to: 

k  =  k  (1  +  y(T  -  T  )]  =  k  (1  +  y6) 
o  o       o 


where 


Y  -  8k/8T  =  A, 


and  9(x,t)  =  T(x,t)  -  T  (x) 


Then,  Equation  (62)  becomes 


k  =  1  +  Ak   +  yf 
o 


(61) 


(62) 


(63) 
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having  neglected  the  terra  yAk  8.   The  multiplication  factor 


k[6(y,T)]  =  1  +  AkQ  +  Y9(y,T) 


can  now  be  written  as: 

(64) 

Y  is  the  temperature  coefficient  of  reactivity,  which  could 
be  either  (+)  or  (-) . 

b.   Feedback  Model  #2 

The  temperature  dependence  of  the  multiplication 
factor  can  be  expressed  in  a  general  form,  adding  the  con- 
tributions due  to  the  Doppler  effect,  due  to  structural  ex- 
pansions and  due  to  other  effects,  as: 

(65) 


dk       ,  1N  n       ,T  Nm  , 

dY  =  Vt"-)    +  Vi^    +  y 


where  D  and  E  stands  for  doppler  and  expansion  respectively. 

k  can  be  expressed  explicitly,  using  the  initial  condition, 

k  =  k   at  T  »  I  ,  thus 
o  o 


k  =  1+Ak  + 


_JL_  [T(-^)n  -  T  (— )nl  + 
o   (1-n)  U^T  ;      o4  ;  J 

o 


i+m     a. 


o  a«        o       o 


(66) 


when  n  =  1,  the  third  term  on  the  right  hand  side  becomes: 


YD  a   •  in  - 


(66a) 


It  is  noted  that  a,  =  a„  =  300°K  and  T  is  usually  expressed 
in  °K. 
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III.   STABILITY  ANALYSIS 

A.   POINT-KINETICS  MODEL 

The  reactor  under  consideration  for  this  analysis  has 
the  following  characteristics: 

(a)  Thermal  reactor. 

(b)  Homogeneous,  bare,  slab  reactor. 

(c)  One-group  delayed  neutron. 

(d)  Constant  power  removal. 

(e)  Stationary  -  full  reactor. 

and  normal  operating  conditions  (no  accidents)  are  assumed. 

The  governing  equations  for  this  system  are: 

£P(t)  -  [(l-B)Ak  -B]  P(t)  +  (1+Ak)3  P(o)       (10) 
C6(t)  =  P(t)  -  (P(t)  (21) 

using  the  following  feedback  models: 

(a)  Ak  =  Ak   +  ye  (25) 

(b)  Ak  -  Ak   +   -^-  [T<^|V  -  T  (^V]      (27) 

o     1-n       1         o   i 

o 

where  Ak   is  a  positive  step  insertion  of  reactivity  by 
external  means,  such  as  control  rod  motion.   The  initial 
conditions  for  the  system  are: 

(a)  P(o)  =  (P(o) 

(b)  9(o)  =  0 

(c)  9(o)  =  0,  T(o)  =  T 

o 

For  the  constant  power  removal,  (p   (t)  becomes  P(o)  and 
Equation  (21)  can  be  written  as 

c6(t)  =  P(t)  -  P(o)  (21a) 
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Assuming  P(o)  =  1  and  introducing  a  d imens ionless  power, 

Equations  (10)  and  (21a)  can  be  expressed  as 

£P(t)  =  [(l-3)Ak  -3]P(t)  +  (1+Ak)3  (10a) 

c6(t)  =  P(t)  -  1  (21b) 

1 .   No  Delayed  Neutrons,  No  Reactivity  Input 

The  point  kinetics  model  when  3  =  0  and  Ak   =  0  is 

o 

reduced  to 

£P  =  yd?  (67) 

c6  =  P  -  1  (68) 

After  some  algebraic  manipulations,  Equation  (67)  can  be 

expressed  as 


L  In  P  =  li 
dt       i 

"Cross  multiplying"  Equations  (68)  and  (69)  yields 
*    d  Yc   ' 

p  -  ^-  in  p  -  t1  ee  =  o 

dt  I 


(69) 


(70) 


which  can  be  written  as 
d 


1  Yc   2 
t(P  -  In  P  -  ijS  9  )  .  0 


(71) 


The  expression  inside  the  parentheses  can  be  called  V,  then 


V  =  P  -  InP  - 
dV 


2  T   9 


(72) 


and  therefore     -1-   =    0. 

d  t 


For  V  to  be  a  Liapunov  function,  it  has  to  fulfill  the  fol- 
lowing conditions: 

(a)  V(P,6)  >  0  for  P>0  and  9>0 

(b)  V(l,o)  =  0  (73) 

t 

(c)  V(P,6)  <  0  for  P>0  and  6>0 
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If  V  of  Equation  (72)  is  modified  as 

V  =  (P  -  In  P-l) 


1  Yc  .2 


2    1 

it  is  seen  that  Equation  (71)  is  not  altered,  and  for  a 
negative  reactivity  input  Equation  (74)  becomes 


(74) 


V  -  (P  -  In  P-l)  +  -| 


Y 


c  „2 


Jt 


(75) 


The  expression  (P  -  In  -  P-l)  >  0  for  all  values  of  P  >  0. 
Then  V  satisfies  condition  (73a)  for  all  values  of  P  and 
6  >  0,  when  y    is  negative;  for  a  positive  value  of  y ,  con- 
dition (73a)  is  satisfied  whenever 

1  Yc   2 

(p  -  in  -  p-i)   >  |  —■  ez 

Condition  (73b)  is  clearly  satisfied  for  the  steady  state 
solution  P  =  1  and  0=0.   Condition  (73c)  gives  V  =  0  for 
both  cases  of  positive  and  negative  y>  which  means  that 
asymptotic  stability  has  been  ruled  out.   A  similar  result 
was  found  by  Ergen  and  Weinberg  (24)  using  the  Hamiltonian 
approach.   Due  to  its  simplicity,  this  case  has  been  added 
to  this  work  to  acquaint  the  reader  with  a  general  view  of 
the  necessary  steps  in  the  Liapunov's  Method  formulation. 
This  case  does  not   represent  any  meaningful  practical 
si tuat ion . 

Typical  data  for  a  thermal  reactor  is  obtained  from 
Solomon  and  Kastenberg  (2_6 )  ,  as  shown  in  Table  1.   From 

this  data,  I    is  found  to  be  1.235  x  10   sec.   An  average 

3 
heat  capacity  for  liquid  light  water  is  used,  1  cal/cm  -  °C 
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TABLE  1 

v  Ef 

0.49780 

-1 

cm 

a 

0.261 

-1 
cm 

V 

2.21  x  105 

-1 

cm-se  c 

6 

0.0064 

- 

D 

0.409718 

cm 

a 

1.41 

- 

A 

0.080 

-1 
sec 
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The  case  of  negative  y  has  been  found  to  be  stable 
for  all  values  of  P  and  0,  but  the  stability  is  only 
neutral  (V  =  0) . 

The  case  of  positive  y  is  represented  in  Figure  6, 
in  which  the  system  is  seen  to  be  unstable  for  most  of  the 
operating  ranges  of  P  and  G. 

2 .   One-Group  Delayed  Neutron,  No  Reactivity  Input 

In  this  case  3  is  included  in  the  analysis  and 
Ak   =  0,  the  governing  equations  become 

£P  =  [(l-3)Ak-3]P  +  (1+Ak)3  (76) 

c6  =  P  -  1  (77) 

Ak  =  yQ  (78) 

Substitution  of  Equation  (78)  in  Equation  (77)  yields 

£P  =  [(1-B)Y6-0]P  +  (l+y0)3  (79) 

"Cross  multiplying"  Equation  (79)  and  Equation  (77)  gives, 
after  some  algebraic  manipulations 


P-  _  In  P 


I 


(1-3)66-  —  6+1-  (1+y6)  -^r^" 


£ 


(80) 


Thus 

^[P-InP-  \   f(l-3)Y62+  jS-   9]  -|  (1+ye)  ^^ 


(81) 


Then 

V  =  (P  -  In  P-l)  -ff1  (l-3)62  +  j± 


(82) 


^  -  f  Cl+ye)  ^ 


(83) 


It  can  be  seen  that  the  inclusion  of  the  delayed  neutrons 
complicates  the  Liapunov  function,  obtained  by  the  same 


A3 


\^/, 


2.0 


l.S- 


1.0- 


0.5- 


W//////A 

Stable 


y|  =  10   /°C 


0.5 


1.0 


> 


6(°C) 


Figure  6.   Stability  domains  for  the  case  of  no  delayed 
neutrons,  no  Ak   input,  positive  temperature  coefficient 

r  O 

of  reactivity. 
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steps  as  for  the  previous  case,  but  this  case  approaches  a 
more  real  situation. 

Two  cases  are  analyzed: 

a.   Negative  reactivity:  y    =    -     \y\ 

B 


v  =  (p  -  in  p-i)  +  |  ^p-L  (i-B)e2 


+  T^e 


(84) 
(85) 


It  is  seen  from  Equation  (84)  that  V  >  0  for  all  positive 

values  of  P  and  G.   Normally  (l-|y|0)  >  0,  then  for  values 

of  0  <  P  <  1,  V  is  less  or  equal  than  zero.   Clearly  V(1.0) 

=  0,  then  V  only  vanishes  at  the  equilibrium  state.   It  is 

concluded  that  asymptotic  stability  is  obtained  whenever  P 

be  less  than  1.0,  and  for  all  9  >  0.   Since  Ak   =  0,  P  can- 

o 

not  be  larger  than  1.0.   Figure  7  shows  the  domain  of 
stability . 

b.   Positive  reactivity:  y    =  |y| 
Equations  (82)  and  (83)  represent  this  case.   The  condition 
of  asymptotic  stability  is  obtained  whenever 

(P  -  In  P-l)  +  j£  9  >  \      |  (1-3)Y92 

and  P  <  1. 

This  case  presents  a  very  narrow  region  of  stability,  as 

seen  in  Figure  8. 

3  .   One-Group  Delayed  Neutron,  Step  Reactivity  Input 

The  two  previous  cases  have  only  theoretical  inter- 
est because  they  do  not  represent  a  practical  situation. 
The  case  of  Ak   input  represents  a  real  practical  situation, 
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p 

1 

■ 

y 

1 

l/l 

III  Ill/Ill  """Ill/ill 

1/1/    1 //Stable  ////// 

II I'll II hu,  01  ll", 

w 

0 

e 

Figure  7.  Stability  domains  for  the  case  of  delayed 
neutrons,  no  Ak  input,  negative  temperature  coeffi- 
cient of  reactivity. 
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Figure  8.  Stability  domains  for  the  case  of  delayed 
neutrons,  no  Ak  input,  positive  temperature  coeffi- 
cient of  reactivity. 
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i?    m    [(1-3) (Ak  +  y0)-3]  P  +  (1+Ak  +  y0) 


i.e.  start-up  and  shut-down,  planned  or  unplanned  motion  of 
control  rods.   The  governing  equations  for  this  system  are: 

(86) 

c6  -  P  -  1  (87) 

Substituting  Equation  (87)  in  Equation  (86)  yields,  after 
simplifying  terms: 


P  -  [Ako(l-g)-3]  j 


.1 


£(AkQ  +  3y6)  +  p-  (l-3)P 


(88) 


Then 


-fp-l-[Ako(l-3)-3]  f  el  =  |(Ako+3yG)H-  J(1-3)6P       (89) 


Thus 


V  =  (P-l)  +  [3-  AkQ(l-3)]  j   ' 
v  =  x<Ako  +  3ye)  +  *  (1-3)6P 


(90) 
(91) 


For  negative  y>  Equation  (91)  remains  unchanged  and  Equation 
(90)  becomes 


V  =  f(AkQ  -3|Y|  e)  -  ■Lfi—  (1-3)6P 


xLa- 


(92) 


In  order  to  obtain  a  positive  definite  V,  it  is  required 
that 


P  +  [3  -Ak  (1-3)]  y  6  >  1 

O  36 


This  condition  is  satisfied  whenever 

8 


Ak   < 


(93) 


(94) 


"o   1-3 

for  all  6  <  0 ,  in  agreement  with  the  linear  theory.   To  get 
a  negative  definite  V,  it  is  required  that 

Ak  <  [3  +  (l-g)P]  |y|9  <  0  (95) 

This  is  the  additional  stability  condition  required  by  the 
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nonlinear  theory.   Equations  (94)  and  (95)  are  specifying  the 
bounds  for  asymptotic  stability  for  a  negative  y.   It  is 
noted  that  immediately  after  the  step  change  in  the  multipli- 
cation factor,  the  neutron  flux  increases  rapidly;  this  is 
referred  to  as  the  prompt  jump,  Lamarsh  (25),  and  it  is  ex- 
pressed as 


P(t) 


3(1-Ak  ) 

r 

"0-2'k" 


P(o) 


(96) 


For  positive  y>  it  i-s  seen  from  Equation  (91)  that  V  will 

always  be  positive  definite  for  all  values  of  P  and  0, 

representing  an  unstable  situation.   Data  from  Table  1  is 

used  to  specify  domains  of  stability  for  various  Ak   inputs. 

The  first  bound  for  all  these  cases  is  provided  by  Equation 

(94)  giving  a  Ak   <  1.00625$,  in  order  to  obtain  a  positive 

definite  V,  for  all  P  >  1  and  0  >  0.   It  is  clearly  seen 

that  P  has  to  be  greater  than  one,  the  steady  state  power, 

after  a  Ak   insertion.   The  second  bound  is  obtained,  set- 
o 

ting  V  =  0  in  Equation  (92),  Figures  9  and  10.   After  the 
insertion,  the  power  will  increase  to  an  amount  specified 
by  Equation  (96).   From  Figures  9  and  10,  it  can  be  seen 
that  right  after  the  Ak   insertion,  the  system  is  unstable, 
and  P  and  0  increase,  until  the  boundary  V=0  is  crossed, 
then  the  system  becomes  asymptotically  stable  and  the  tra- 
jectory returns  to  the  equilibrium  position,  but  when  V=0 
is  crossed  again,  the  system  turns  to  be  unstable,  and  in 
that  way,  the  trajectory  keeps  oscillating  along  the  line 

• 

V=0.   A  similar  situation  is  obtained  by  Meghreblian  and 
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Ak   =  0.156$ 
o 

y   =  10~4/oC 


|  Prompt  jump 


10 


100 


— j 

1000 


8(°C) 


Figure  9.   Stability  domains  for  the  case  of  delayed 
neutrons,  step  Ak   =  0.156$,  negative  temperature 
coefficient  of  reactivity. 
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1.45,. 


Ak  =  0.3125$ 

?     -4 
Y  =  10   /°C 


Prompt  Jump 


1000 


6(°C) 


Figure  10.   Stability  domains  for  the  case  of  delayed 
neutrons,  step  Ak   =  0.3125$,  negative  temperature 
coefficient  of  reactivity. 
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Holmes  (2_2_)  when  solving  the  system  of  differential  equa- 
tions.  Their  results  showed  that  this  is  the  case  of  damped 
oscillation,  the  presence  of  the  delayed  neutrons  providing 
the  damping  force,  thus,  following  an  insertion  of  reacti- 
vity, the  system  will  reach  a  new  steady  state  power. 

4 .   One-Group  Delayed  Neutron,  Ramp  Reactivity  Input 

The  governing  equations  for  the  system  are  Equations 
(86)  and  (87),  with  the  difference  that  now  Ak   is  a  func- 


tion of  t  expressed  as 


Ak  (t)  =  at 
o 


(97) 


where  a  represents  the  rate  of  insertion  of  reactivity. 
Combining  Equations  (86)  and  (87)  yields,  after  some  alge- 
braic manipulations 

li    -((l-3)yc  60+  gee  -<l-3)atP  -  3at  =  ye  (98) 

Then 

d  1        Yc   2    ec 

^[p   -i-l  (i-6)   ^e2   +-f  6]  - 
=  (1-3)  ~  P  +  p  +  3  f^- 
From  Equation  (99),  V  and  V  are  obtained  to  be: 

v  =  p  -  i  -  |  (i-3)  p  e2  +  j±  e 


(99) 


(100) 


V  -  (1-3)  4  P  +  f  +  fip 
For  a  negative  y>  Equations  (100)  and  (101)  become 


(101) 


V  =  P  -  1  + 


3c 


yc6 


+  f  (1-6)  t 


at 


V  =   i_b  z±.   p  _  X6-  + 


at 


(102) 
(103) 
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The  following  condition  is  obtained  for  asymptotic 
stability : 

Ak  (t)  =  at  < 


_xl 


(104) 


(l-3)P+3 

It  can  be  seen  that  V  is  positive  definite  for  all  P  >  1 
and  8  >  0.   Equation  (104)  is  necessary  to  obtain  a  negative 
definite  V  to  ensure  asymptotic  stability.   For  a  positive 
Y,  the  system  is  inherently  unstable.   Equation  (101)  gives 
a  positive  definite  V  for  all  values  of  P  and  6  for  t  >_  0. 
A  typical  ramp  insertion  of  reactivity  is  shown  in  Figure  A. 
During  the  reactivity  insertion,  the  system  is  unstable, 
then  V  =  0  is  obtained  after  the  ramp  insertion  has  ended 
(10  sec)  and  the  power  and  temperature  have  risen  to  an 
appropriate  level.   Setting  V  =  0,  in  Equation  (103),  the 
following  expression  is  obtained 

|  -y  |  S  —  3Ak 


P  = 


(l-B)Ak 


(103a) 


At  the  end  of  the  Ak   insertion  (t  =  10  sec),  this  expres- 
sion implies  that  P  depends  linearly  on  6,  Figure  11. 

5  .   One-Group  Delayed  Neutron,  General  Reactivity  Input 
The  equations  for  this  system  are: 

&P  =  [  (l-3)Ak-B]P  +  (1+Ak)3  (106) 

c6  -  P  -  1  (107) 

Y, 


Let 


o 


R<T>-  3^  [T(^V  -  To(^°)n] 

o 


(108) 


(109) 
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Ak      =    O.OOOlt 
o 


Figure    A 
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Ak  =0.156$(10  sec) 

?     -4 

Y  =  10  7°c 


1000 


Figure  11.   Stability  domain  at  the  end  of  a  ramp 
insertion  of  reactivity. 
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Then  Equation  (106)  becomes 

£P  =  [(l-3)(Ak  +Yn-R(T)-3]P  +  [1+y  -r(t)]B     (110) 

o   D  D 

Combining  Equations  (110)  and  (107)  yields,  after  some 
algebraic  manipulations 


I?    -  (l-3)Ak  c0+3c6  -  [  (1-3)P  +  3]YT,R+Ak 


D 


(HI) 


Then 


^tp-1  + 


3-(l-3)Ak 

_ 


c0]  = 


Yn'R    Ak 
=  [(1-3)P  +3]  ~—  +  -~ 


Therefore 


v  -  p-l  +  [3  -  (1-3)  Ak  ] 

o 


c9 


and 


YD'R 
V  =  [(1-3)P  +3] 


Ak 


(112) 


(113) 


(114) 


I  "   I 

Equation  (113)  is  independent  of  y  ;     therefore,  V  is  positive 
definite  whenever  conditions  (93)  and  (94)  are  satisfied. 
The  similarity  with  the  step  insertion  of  reactivity  can  be 
seen.   Equation  (114)  becomes,  for  negative  y 
Ak 


V  = 


IydI-r(t) 

-  [3  +  (1-3)P]  - 


Equation  (115)  is  negative  definite  whenever 
AkQ  <  [3  +  (1-3)P]  |ydI 'R(T) 


(115) 


(116) 


for  all  n  ^  1.   The  case  n  =  1  is  the  logarithmic  case  and 
Equation  (116)  becomes 


Ak   <  300[3+(l-3)P]   Yj  In  ^- 
o  D       T 

o 


(117) 
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The  case  of  positive  y    leads  to  an  unstable  situation  due 
to  the  fact  that  Equation  (114)  is  positive  definite  for 
all  n,  p  and  T.   Data  for  this  case  are  found  in  Thompson 
and  Beckerley  (2J3)  and  are  typical  of  a  fast  reactor.   The 
following  specific  case  is  analyzed: 

(a)  Oxide  reactor  with  volume  ratio  U0o  to  P  0.    =    7 

I  u  2 

(b)  Sod j urn  density  =  50% 

(c)  YD  =  1.28  x  10~5/°K 

(d)  n  -  0.96 

(e)  3  -  0.0033 

(f)  Ak  -  0.15  6$ 

o 

The  domain  of  stability  is  shown  in  Figure  12,  where  the 
V  =  0  is  obtained  from  Equation  (115). 

B.   DISTRIBUTED-PARAMETER  REACTOR  SYSTEM 

The  reactor  under  consideration  for  this  analysis  has 
the  following  characteristics: 

(a)  fast  reactor 

(b)  homogeneous,  bare,  slab  reactor 

(c)  one-group  delayed  neutron 

(d)  Newton's  law  of  cooling 

(e)  stationary-fuel  reactor 

and  normal  operating  conditions  (no  accidents)  are  assumed 
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1200 


e(°K) 


340 


Figure  12.   Stability  domains  for  the  case  of  delayed 
neutrons  and  general  reactivity  insertion. 
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The  governing  equations  are: 

2 


3  *(y»T)  +[Cl-B)ka-l]»Cy.T)  +  jg  X^(y,x) 

9y2  VZf 


„  Wy,-Q 
vLa  — ar     =  -X-f  (y,x)  + 


3t 


Pg 


»P(y»T) 


(37) 


(46) 


c,vE   86^>T)  =  E  Ef  ^(y,T)  -  h6(y,T) 
r   a    oT  r 


(57) 


using  the  following  feedback  models: 

(a)  k(y,T)  =1  +  AkQ  +  Y6(y,T) 

(b)  k  =  l+Ako+  Y^-[T(^)n  -  Vf^"1  + 

o 

+  ^-[T(^-)m  -  T  (^)m]  +  Y  (T  -  T  ) 
1+m    a~       o  a„        o       o 


(63) 


(66) 


where  Ak   is  an  external  positive  reactivity  insertion.   The 
o  r  J 

boundary  conditions  for  the  system  are: 

(a)   i|>(+  W,  t)  =  0 
The  initial  conditions  are: 

(a)  iMy,o)  -  0 

(b)  £(y,o)  =  0 

(c)  0(y>o)  -  0 

The  case  n  =  1  will  give  the  logarithmic  term  in  the  doppler 
expression . 

1 .   One-Group  Delayed  Neutron,  Step  Reactivity  Input 

To  formulate  the  concept  of  stability  a  distance  be- 
tween the  equilibrium  position  (d>  ,  C  ,  T  )  and  any  perturbed 

ooo 
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position  ((J),  C,  T)  is  defined.   The  variation  of  this  norm, 
which  is  induced  by  the  inner  product  of  the  Hilbert  space 
in  which  the  solutions  of  the  system  are  defined,  with  time 
will  provide  information  concerning  the  perturbed  position. 
Let  W 


[d(4>,(}>o;C,Co;T,To)]2  =  J      U±2*2+    A^  2+A326 2)  dy 


where 


and 


-W 


e2Zf2 


e:2v2Z  2 


232v2 


*  2    „  2  2„  2 
A„   =  Cr  v  E 
3      fa 


\\>    =  <f>  -  <J)o 

•d    =  C  -  C 

t  o 


(118) 


(119) 
(120) 

(121) 


The  Liapunov  function  is  then  defined  as 


v(<k  g,  e)  = 


=  <  d  ,  d  > 


(122) 


which  by  definition  is  positive  definite.   It  is  noticed 
that  the  introduction  of  the  constants  A..,  A?  and  A„  makes 
V  represent  the  total  rate  of  energy  increase  per  unit 
volume  of  the  system,  which  is  defined  as  the  distance  be- 
tween the  equilibrium  and  the  perturbed  states. 

V  vanishes  only  at  the  equilibrium  position  (<J>  ,C  , 
T  );  therefore,  the  requirements  of  a  negative  definite  V 
will  provide  the  domain  of  asymptotic  stability. 
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Thus 


W 
di         /    ^    Lf    ^3y    +         o2„2    *      3y    "     *""*     v    "a         8y 


Ce    Z.    ip-^1   +    a    <£     ~   +    2C^    v    I       6    -r— )     dy 

2      ?  "v  fl  hv 


-W  (123) 

Substituting  Equations  (37),  (46)  and  (57)  in  Equation  (123) 
yields 


W 


g-    E2Zf2 
dT  f 


■W 


#  ^2   +    [(l-B)ka-l]^2   +  ^  £  ip 


3y 


vZ      ka 

+  __£ —  4$  _ 

3v2Zf2 


vE    A  -         2C,Z      ka  2C£vZ    h      „ 

a  /2.  fa  „,  f       a       ~  2        . 

£     +  : —  ^  "  — : — : —  e      dyi 


32v2zf2 


veE 


£2E    2 


(124) 


Integration  by  parts  of  the  first  term  and  using  the 

boundary  condition  i|K+  W,  t)  =  0  yields 
W  W 


(125) 


-W  -W 

Knops  and  Wilkes  (_27_)  ,  in  their  work,  provide  the  inequali- 
ties useful  for  this  kind  of  analysis.   The  so-called 
"eigenvalue  inequality"  plays  an  important  role  in  this 
specific  problem.   For  a  system  defined  by  the  eigenvalue 
problem 

V2u  +  Ay  =  0 
in  domain  with  u(a,t)  =  y(-a,t)  =  0,  on  the  boundary,  the 
following  inequality  can  be  stated: 


liZ  dy 


*/ 


^>2  * 


(126) 


-a 


-a 
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an  inequality  that  can  be  proved  using  the  calculus  of 
variations.   (See  Appendix  A).   In  applying  this  inequality 
to  this  work,  it  is  assumed  that  the  perturbed  reactor  has 
the  fundamental  eigenvalue  not  appreciably  different  from 
that  of  the  unperturbed  reactor,  then 


4W' 


(127) 


Thus,  Equation  (126)  becomes 


w 


V  dy 


(128) 


-W  -W 

Substitution  of  this  result  in  Equation  (124)  yields 


dt  — 


W  / 


e2Z 


B2L2ip2  -  [(l-$)ka-l]Tj/2 


.  ,        vi   ka 


f  * 


+  —2 f>2 

a2v2zf2    c 


2n  v£   ka         2CrvZ  h 

f  a e*  + 


veZ 


e2£  2 


dy 


(129) 
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After  introduction  of  the  feedback  term,  the  following 
expression  is  finally  obtained: 


dV  .  2   v   2 

dx  —        f 


[l-(l-3)(l+Ako)]  f 


-W  VZ 

+  (1-3)Y6^2  -  (l+Ak  )(-£-+ — 5 

0   Vlf    3v2Zf2 


,     vE 
+  (4-  + 


■vZ 


3v2Zf2 


)  yegip 


vZ  X 


S  v2Z  2a 
f 


^2    * 


2C  vE  h       2C  vZ 
+  e2 1 §_  (1+Ak  j  0^  + 


e2Z  2a 


cZ.v 


2C  vE 

+  a   Y62^   I  dy 

eZ  2v 


(130) 


Let 


a-  -  1  -  (1-3) (1+Ak  ) 
1  o 


vZ  X 
a 


*2    32v2Zf2a 

2CrvZ  h 

f   a 
a   = 

3    e2Zf2a 

1     vE  ! 
X      a_ 


4   vEf    3v2Zf2 


a5  = 


2C.vZ  ' 
f   a 

eZf  2v 


(131) 
(132) 

(133) 
(134) 
(135) 


62 


Then  Equation  (130)  can  be  expressed  as 

W 


4^  <  -  e2  E.,a  a 

dx  —        f 


2  +  a3  e2   + 


/  [ax  ^2  +  a2^ 

-W 
+  (1-3)  yQ^2  -  (1+Ak  )  a4^>  +  a-4  Y©£*  - 


-     (1+Ak  )  a5  B\p    +  a   Y^ll  dy 


(136) 


To  obtain  asymptotic  stability,  it  is  necessary  to  prove 
that  the  expression  inside  the  curly  brackets  is  positive 
definite.   A  strategy  successfully  employed  by  Buis  and 
Vogt  (2_8)  will  be  used  here.   Equation  (136)  can  be  written 
as  f ollows : 


dV 
dx 


!W 
[      J     (a±    i>2    +    a2£2    +    a3    02)dy]     ■ 


-W 


w 


1- 


J  [  (l+Ako)aA£ij;-a4Ye^-(l-3)Ye^2  +  (l+Ak    )  a    Gip-a    y9  2iJj]  dy 


-W 


r w 

J      (ax    i>2    +    a2    ^2    +    a3    62)dy 


(137) 


(138) 


-W 

Equation  (137)  can  be  expressed  concisely  as 

4^-  <  -  R  •  P  *  Q 

dx  —  x 

where  R  =  e  Z   a  is  essentially  a  positive  constant, 


fW  2 

P    =      /         (a       i>2    +    a2  ^       +    a~    02)     is    greater    than    zero 

'  -W 
whenever  a..  a~  and  a»  are  greater  than  zero.   It  is  already 

known  that  a~  and  a„  are  positive  constants. 


Thus 

zero  in  order  to  have  P  >  0. 


a-  =  1  -  (1-3) (1+Ak  )  has  to  be  greater  than 
1  o 
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From  this  condition,  the  first  bound  for  stability  is 
ob  t ained 


Ak      < 

3 

o 

1-3 

(139) 


in  agreement  with  the  linear  theory.   To  prove  that  Q, 

second  expression  in  brackets  in  Equation  (137),  is  positive 

definite,  the  Buniakowsky-Schwarz  inequality  will  repeatedly 

be  used 

a  a  a 

fig    dy  <  (  /  |f  |2dy)1/2(/|g|2dy)1/2    (140) 
-a  -a  -a 

Due  to  the  positive  Ak   insertion,  the  system's  state 
r  o 

variables  ty , £ ,  6  are  positive  for  any  t  >  0;  therefore, 
the  absolute  value  in  the  inequality  is  not  necessary.   It 
is  also  known  that 


W 


2     ,  2 


•W 


f    A±       ty       dy    < 


W 


2^2 


_W/Vf    d 


13 12 


w 


-w 


/ 


2    .2 


dy    < 


(141) 
(142) 
(143) 


State  variables  in  engineering  have  the  properties  of  func- 
tions in  the  Hilbert  space,  and  the  norm  in  this  space  is 
specifically  defined  as 


(  /  y2dx) 


1/2 


=    y 


(144) 


-a 


Due  to  the  fact  that  the  state  variables  in  the  present 
work,  \p ,  ■£>  and  9,  are  continuous  in  space  and  time  domains 
together  with  their  derivatives  to  arbitrary  kth  order,  it 
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is  expected  that  they  also  have  the  properties  of  functions 
in  the  Banach  space.   The  norm  in  the  Banach  space  (L  )  can 


be  de f  ined  as : 


(yV  dx)1/m  -  |p|,   for  m    1 


(144a) 


m 


-a 


Then,  repeated  applications  of  these  inequalities  yields 


°*** &+***-*&&:*  ft?  +ih*tf 


1    - 


An        VA„  A    '•'"""  A,  NAnA0  AnA»  .     2 

1  2  3  123  13  A„ 


min( — j     , 

a: 


2  '  2 

A2         A3 


)IU 


>    0 
(145) 


where  the  minimum  is  taken  in  order  to  assure  the  inequality. 
Therefore,  for  Q  >  0,  it  is  necessary  that 


(1+Ak  )  a.  ac 
o  (_4  _5} 

A      A  A 

A!    A2  A3 


al    a2    a3 
min  ( 2>  — 2  ,  2) 

Al    A2    A3 


2Wy  r_h 

Al   A2A3 


A1A3    a/ 


(146) 


This  is  the  additional  stability  condition  prescribed  by  the 
nonlinear  theory.   The  inequality  (146)  gives  the  second 
bound  for  asymptotic  stability  of  the  system,  which  repre- 
sents a  spherical  surface  in  the  first  octant,  i.e.  \p ,  "Jf 
and  9^0.   Domains  of  stability  are  represented  in  Figure 
13.   The  solutions  for  ^(y,T)  and  6(y,x)  can  be  obtained 
from  Equations  (46)  and  (.57),  as  follows: 


J> 


£  (y,T)  -   Jh^t-t')  i>(y,T')  dT1 


(147) 
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Traj  ectory 


W 


>    (A     2     I      rf,*^1'2 


(A 


iT'dy) 


-W 


Figure  13.   Schematic  of  the  spherical  surface  determin- 
ing stability  domains  for  the  distributed  parameter 
reactor  system  after  a  Ak   insertion. 
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where 


A  (t-t') 

vZ 


H  (t-t')  =  vEf3  e 


(148) 


and 


e(y,x) 


H2(t-t')  ^(y,x')  di' 


(149) 


where 


£l 


h  (t-t') 
CvZ 


H2(t-t')  - 


CvL' 


(150) 


Application  of  the  mean  value  theorem,  in  the  time  domain, 

yields 

fT 

(151) 


-£(y,T)  -  iKy,r)  /  H1(T-T')dT'  =  ^(y,x)T1(T) 


where 


BwE  Z        -  -—  t 
TX(T)  =  f-^    (1  -  e  wLa         ) 


(152) 


and  0  <  T  <  T . 


Also : 


(y,T)  =  ^(y,T)  /  H2(t-t')  dx'=  ^(y,T)T2(T)     (153) 


where 


eZ 


T9(T)     -1  (l  -  e   CvZa    ) 


(154) 


and  0  <  t  <  T  . 

Substituting  Equations  (151)  and  (153)  in  Equation  (118) 


results  in 

II  Til  2 


A12^2(y,T)dy  +  /  A22^2(y,T)T12(T)dy  + 

w 


A32^2(y,T)  T22(t)  dy 


(155) 
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Applying  the  mean  value  theorem,  in  the  space  domain,  yields 


lid 


Ax2iJ;2(y,T)  +  2WA2V(y,T)  T^CO 


+  2WA3V(y,T)  T22(x) 


(156) 


or 


d||      =  i|;(y,T)    v/a12+  2wa2?t12(t)+2wa32t22(t) 


(157) 


where  0  <  1     <    T  and  -W  <  y  <  W. 

Then  the  mean  value  flux  can  be  expressed  explicitly  as 

follows : 


i|>(y,T)  = 


11 


(158) 


r 

•  A. 


+  2WA,  T,(t)  +  2WA,  T0z(t) 


1   '  """2  U  **'  '  """3   2 
Typical  values  for  the  fast  reactor  nuclear  parameters  are 
listed  in  Table  2,  after  Solomon  and  Kastenberg  (2_6)  .   Using 
these  parameters,  the  constants  for  the  system  are  found  to 

be: 

13    2     -  2 
a„  =  3.292  x  10    cm  -sec 

27    -2    -1     -1 
a   =  6.185  x  10    cm   -°C   -sec 

9         -1 
a.  =  9.050  x  10    cm-sec 

a5  =  1.470-  x  1016  cm"2-°C~1-sec~1 

A-  =  2.932  x  10"26cal2-cm~2 

A„  =  2.780  x  10~8  cal2-sec~2 

A32=  1.920  x  102   cal2-cm~6-°C~1-sec"2 

W  -  3.0 


and 


Tl(T)=  70.05  d-e-O.OOOCm^ 
T9(t)=  2.38  x  iO-12Cl-e-°-00004A8T) 


(159) 
(160) 
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From  Equation  (139),  the  first  bound  is  given  by: 


Ak   <  1.00331$ 
o 

From  Equation  (146),  the  norm,  second  bound,  is  obtained 

to  be : 

.5 


||d]|   >   2  x  10' 


for  a  Ak   =  0.156$.   The  mean  value  flux  from  Equation  (158) 

O  IX* 

is  shown  in  Figure  14.   The  physical  meaning  of  this  mean 
value  flux  ip(y,l)  is  depicted  in  Figure  15.   The  mean  value 
flux  provides  an  order  of  magnitude  information  of  the  in- 
crease in  neutron  population  at  the  time  the  surface  of 
stability  is  reached.   From  Figure  14,  it  is  seen  that  a 
fast  rising  flux  will  quickly  reach  the  surface,  whereas  a 
more  slowly-rising  flux  will  take  longer  to  reach  the  sta- 
bility domain. 

2  .   One-Group  Delayed  Neutron,  Space  Dependent 
Step  Reactivity  Input 

This  system  is  a  special  case  of  the  previously 

analyzed  system.   It  is  desired  to  find  an  answer  to  the 

following  question:   Given  a  positive  reactivity  insertion 

Ak  ,  is  it  safer  to  insert  this  reactivity  uniformly  across 

the  reactor,  Figure  16a,  or  to  insert  less  reactivity  in  the 

central  region,  Figure  16b,  provided  the  total  reactivity 

insertion  is  the  same? 


Equation  (129)  can  be  rewritten  as: 
W 


dV 

dx 


<  -e2Ef2a  f    |  [l-(l-3)k]ijj2  + 


vE  X 
a 


2, ,2 


4:1    + 


S^v'E  *a 


2CrvE  h 
f   a   62 
e2Ef2a 


vE 


vE 


(X+ 


BvEf  Y: 


2C  vE   k     ) 
^  i> ^—  6*  }dy 


(161) 


veE 
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Figure  14.   Mean  value  flux  vs.  dimens ionles s  time  plot 
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Point (y,T) 
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Figure  15. 
reactor . 


Mean  value  flux  schematic  in  the  slab 
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Figure  16a 
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Figure.  16b 


Figure  16.   Space-dependent  Ak   insertion. 
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Equation  (161)  can  be  rewritten  as: 

-2W/3        2W/3  U 

|^  <  -e2Zfa  [  /  {-}  dy  +  J  "{-}  dy  +  /    {-}  dy 

-W  -2W/3         2W/3 

Three  steps  Ak   are  introduced  in  the  following  way 


(162) 


k  =  1  +  Ak   -  yi 


k  =  1  +  Ak. 


yO 


k  =  1  +  Ak   -  yi 


for 
for 
for 


2W 
•  W  <  y  <  -  ~ 

2W  „       2W 

—  1  y  1    3~ 

f-   <  y  <  W 


(163) 
(164) 
(165) 


Substituting  Equations  (163),  (164)  and  (165)  in  Equation 
(163)  yields,  after  some  algebraic  manipulations  and  intro- 
duction of  the  constants  a,,  a~,  a~,  a,  and  a,: 


:w/3 


dV  ^     2v  2 
dT  —       f 


[a  %2  +  a2^2  +  a3G2+  (l-3)y9ij): 


-W 


-(1+Ak  )a  ,4^+    a^yQ{:^~    (1+Ak  )a56i|>  +  a  y9  \p]dy+ 


2W/3 


+    j       [a1'>2+    a2^2+    a3e2+(l-3)yei|;2-(l+Ak2)aA^>    + 
-2W/3 

+    auyQ-^-\p    -    (1+Ak2)a3ei|j    +    a5y92^]    dy    + 


/   [ai 


+ 
2W/3 


%2+    a2g  2+    a3    92+    (l-p)y9^2-(l+Ak1)aA^  r\)   + 


+    a4ye^-ijj    -(1+Ak    )a5    Qi>   +    a5    y92ip]dy 
where 


a^      =    1    -    (1-8) (1+Ak1) 


a^'    =    1    -    (1-3)(1+Ak2) 


(166) 

(167) 
(168) 
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Following  the  same  steps  as  in  the  previous  analysis,  the 
following  conditions  for  asymptotic  stability  are  obtained: 


Akl  <  T^ 


(169a) 


Ak2  <  iq 


conditions  which  agree  with  the  linear  theory,  and 
(1+Ak.)  a 


dl|x> 


1'  -4    a5,     .    *1    a2    a3 

a7~(a7  +  a7}   ~  min   (~T'  77'   7T 

12  3  A,  A_  A„ 


) 


w  X_(li 


3    A2    A£A3 


,     (1-3)     ,     "5 

A    A  2} 

A1A3         A. 


(169b) 


(170) 


II a 1! 


'2 


(1+Ak2)/a4         a5,  .        ,al         a2         a3 

~aT~  caT  +  a^   "  mln   (7T'  7T'  72 

12  3  A..  A„  A~ 


pi  r_(!i_  +  (141  +  !_!_) 


3       Ax    A2A3 


A1A3 


(171) 


These  are  the  additional  condition  of  stability  required  by 
the  nonlinear  theory.   In  order  for  the  total  reactivity 
insertion  to  remain  the  same  for  both  strategies,  one  must 


have 


-W 


-2W/3 


2W/3 


ik  J  Akody  ■  I J  Akidy  +  h  J   Ak2d? 

W  ~W  -2W/3 


(172) 


From  Figure  17,  it  is  seen  that  in  Region  2  in  which  Ak^  < 

Ak  ,  (I  d  |L  >  !|d|j,  so  that  the  region  of  stability  is  attained 

sooner  than  in  the  case  of  uniform  Ak   insertion;  but  ||d|L> 

||d||,  so  in  Region  1  of  the  reactor,  where  Ak..  >  Ak  ,  the 

region  of  stability  is  attained  later  than  for  the  case  of 

uniform  Ak   insertion.   However,  the  flux  in  the  central 
o 
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Figure  17.   Spherical  surfaces  for  the  case  of  space' 

dependent  Ak   insertion. 
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region  of  the  reactor  plays  a  more  important  role  in  the 

transient  behavior  of  the  reactor;  consequently,  the  space 

dependent  reactivity  insertion,  with  Ak„  <  Ak  ,  is  safer 

than  the  uniform  Ak   insertion. 

o 

3 .   One-Group  Delayed  Neutron,  General  Reactivity  Input 

For  this  case,  Equations  (37),  (46),  (57)  and  (66) 

are  the  governing  equations.   The  feedback  coefficients,  y> 

Y^  and  Y„ .  are  assumed  to  be  negatives;  Y_  and  Y_  are  the 
D        E  D        E 

Doppler  and  Expansion  coefficients  of  reactivity,  respective- 
ly.  The  same  Liapunov  function,  Equation  (122),  is  used 
and  a  similar  development  previously  discussed  leads  to  the 
following  equation,  in  analogy  to  Equation  (136): 


dV 


■W 


dT    <    -e2Zf2cJ   /  [a1^2+    a2p    +    a302+    (l-$)ij> 2F(0)    + 

+    a4       lf>F(6)     a56^F(6)     -     (l+AkQ)  a^p  i>~  (l+AkQ)  a^ip  ]  dy  I  (  173  ) 


where  the  feedback  term  is 

n 


YtJR 


,m 


F(9)  =  |y|6  +  -T^ [ 


,1-n 


T  i-11]-^  '^-[l1^1-!  1+m](174) 
o  1+m  o 


here  R  and  P  are  constants  equal  to  300°K.   F(G)  is  positive 
for  all  values  of  n  >  0  and  m  >  0.   It  is  noted  that  the 
Doppler  feedback  is  the  primary  feedback  mechanism  in  large, 
ceramic  fueled  fast  reactors,  and  expansion  is  usually  the 
dominant  feedback  in  small  metal  fueled  fast  reactors.   For 
n=l,  Equation  (66a)  is  used.   A  different  approach,  more  con- 
servative, is  used  to  analyze  Equation  (173).   The  expression 
inside  the  curly  brackets  in  Equation  (173)  has  to  be  posi- 
tive definite  in  order  to  ensure  asymptotic  stability.   With 
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this  in  mind,  a  comparison  of  the  orders  of  magnitude  among 
the  terms  is  carried  out. 


dx  —       f 


a{XUi 


^2+${a7$  -    (1+Ak  )a,ij>}  + 


o'"4 


+  6  {a39  -  (l+Ako)a5ip)  +  i];F(e){(l-3)ip+a  f+  a56}]dy|  (175) 

In  order  to  get  a  positive  definite  function  inside  the 
curly  brackets,  the  following  conditions  have  to  be  satis- 
fied: 

o 

(a)  a..  >  0  which  yields  Ak   <  -~e 

1  J  o    1-p 

(b)  a2'f^  -  (1+Ak  )  a4^  >  0  which  yields 


and 


i  <  zi 

f>  (l+Ako)aA 


(c)       ao0    -     (1+Ak    )ac^    >    0    which    yields 
3  o      5  J 


ik 


6  (1+Ak    )ac 

o        5 

Typical  values  for  the  feedback  coefficients  are: 

Y=  10_A/°K,  v  =  7xlO"6/°K,  y    =    9.5  x  10~6/°K.   Evaluating 
1)  b 

these  conditions  for  a  typical  fast  reactor,  with  a  Ak   = 
0.156$,  it  is  obtained  that: 

(a)  Ak   <  1.00331$ 

o 

(b)  t   <    3634.0 

A  ii 

(c)  I  <  4.203  x  1011 

These  three  conditions  have  to  be  satisfied  simultaneously. 
The  domain  of  stability  is  shown  in  Figures  18a,  18b  and 
18c.   It  can  be  seen  that  the  domain  of  stability  is  very 
large  when  the  coefficients  of  reactivity  are  negative. 
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Figure  18a.   Stability  domains  for  the  distributed 
parameter  reactor  system  after  a  general  reactivity 
ins  er t ion . 
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Figure  18b.   Stability  domains  for  the  distributed 
parameter  reactor  system  after  a  general  reactivity 
insert  ion . 
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Figure  18c.   Stability  domains  for  the  distributed 
parameter  reactor  system  after  a  general  reactivity 
ins  er t ion . 
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IV.   CONCLUSIONS 

Liapunov's  Second  Method  has  proven  to  be  a  useful 
method  to  find  stability  domains  in  Nuclear  Reactor  Control. 
The  fact  that  a  Liapunov  Function  is  not  unique  makes  the 
choice  of  such  a  function  to  depend  upon  the  experience  and 
ingenuity  of  the  analyst.   Once  an  appropriate  choice  of  the 
Liapunov  Function  has  been  made  practical  results  can  be 
found  . 

For  the  Point-Kinetics  Model,  a  Liapunov  Function  vas 
obtained  without  imposing  an  a  priori  bound,  as  most  of 
the  conventional  methods  do.   Domains  of  stability  are  pre- 
sented for  different  cases  using  typical  data  from  thermal 
reactors . 

The  Dis t r ibut ed-Parame ter  Reactor  System  was  studied 

using  the  concept  of  a  norm,  an  a  priori  bound,  and  the 

variation  with  time  of  this  norm,  defined  in  a  Hilbert 

space,  was  analyzed  to  obtain  a  domain  of  stability.   Data 

from  a  typical  fast  reactor  was  used  to  define  the  spherical 

surface  of  stability  and  to  estimate  a  mean  value  flux,  in 

space  and  time,  when  this  surface  is  reached.   The  case  of 

space  dependent  Ak   insertion,  with  less  reactivity  in  the 

central  region,  was  found  to  be  safer  than  the  case  of 

uniform  Ak   input, 
o 

The  linear  theory  stability  condition  appeared  through- 
out the  analyses  giving  confidence  to  the  obtained  results, 
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in  addition  to  that  condition  and  additional  condition(s) 
appeared  in  each  case,  required  by  the  nonlinear  theory. 

The  inclusion  of  the  multigroup  delayed  neutrons  is 
recommended  for  further  improvement  of  this  analysis. 


82 


APPENDIX  A 
THE  EIGENVALUE  INEQUALITY 


Let  A  >  0  be  the  least  eigenvalue  of  the  system 

i2. 


4-4  +  X]p  m   0 

dy2 
with  boundary  conditions:   (a)  tJj(+W)  =  0 
the  following  inequality  can  be  stated: 


(Al) 


X     J    i|;2dy   <    /  (^)2  dy  (A2) 

Proof  of  this  inequality  can  be  done  using  tbe  ca 1  cuius 

of  variations . 

Let 


W 

J  (ip'  2-  A^2)dy  >  0 


~W 


Applying  Euler  s  equation:   -: —  F  ,=  F 


dy   ijj'    \p 
Equation  (Al)  is  obtained. 
Solving  the  eigenvalue  problem  A   is  obtained  to  be: 


(A3) 
(A4) 


Then 


A  = 


w 


2.  2 


=  B  L 


(A5) 


4W 


/• 


w 


,dU)x  2  ,      2  2 
<d7>  dy  >  B  L 


^2dy 


(A6) 


■w 
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